#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Generate Merged Map                                                #
#                                                                           #
# 2dx.org, GNU Plublic License.                                             #
#                                                                           #
# Created..........: 02/20/2006                                             #
# Last Modification: 09/20/2006                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 85
#
# MANUAL: This script makes use of MRC and CCP4 commands to generate the final projection map.
#
# PUBLICATION: 3D reconstruction of two-dimensional crystals: <A HREF="http://www.ncbi.nlm.nih.gov/pubmed/26093179">Arch Biochem Biophys 581, 68-77 (2015)</A>
# PUBLICATION: 3D Reconstruction from 2D crystal image and diffraction data: <A HREF="http://dx.doi.org/10.1016/S0076-6879(10)82004-X">Methods Enzymol. 482, Chapter 4, 101-129 (2010)</A>
# PUBLICATION: 2dx - Automated 3D structure reconstruction from 2D crystal data: <A HREF="http://journals.cambridge.org/action/displayAbstract?aid=1943200">Microscopy and Microanalysis 14(Suppl. 2), 1290-1291 (2008)</A>
# PUBLICATION: 2dx_merge - Data management and merging for 2D crystal images: <A HREF="http://dx.doi.org/10.1016/j.jsb.2007.09.011">J. Struct. Biol. 160(3), 375-384 (2007)</A>
#
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: SYM
# DISPLAY: ALAT
# DISPLAY: tempfac
# DISPLAY: avrgamphsNUMBER
# DISPLAY: avrgamphsRESOL
# DISPLAY: realang
# DISPLAY: realcell
# DISPLAY: diffmap_doit
# DISPLAY: diffmap_source
# DISPLAY: diffmap_weight
# DISPLAY: npo_line1
# DISPLAY: npo_line2
# DISPLAY: calculate_subvolume
# DISPLAY: mask_subvolume_PDB
# DISPLAY: mask_subvolume_PDB_file
# DISPLAY: mask_subvolume_PDB_radius
# DISPLAY: merge_alsoevenodd
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set tempkeep = ""
set RESMIN = ""
set RESMAX = ""
set CS = ""
set KV = ""
set ALAT = ""
set realang = ""
set realcell = ""
set rot180 = ""
set revhk = ""
set rot90 = ""
set beamtilt = ""
set tempfac = ""
set SYM = ""
set avrgamphsNUMBER = ""
set avrgamphsRESOL = ""
set det_tilt = ""
set diffmap_doit = ""
set diffmap_source = ""
set diffmap_weight = ""
set merge_ref_num = ""
set merge_comment_1 = ""
set merge_comment_2 = ""
set merge_comment_3 = ""
set merge_comment_4 = ""
set merge_comment_5 = ""
set merge_comment_6 = ""
set merge_comment_7 = ""
set merge_comment_8 = ""
set merge_comment_9 = ""
set merge_register_date_1 = ""
set merge_register_date_2 = ""
set merge_register_date_3 = ""
set merge_register_date_4 = ""
set merge_register_date_5 = ""
set merge_register_date_6 = ""
set merge_register_date_7 = ""
set merge_register_date_8 = ""
set merge_register_date_9 = ""
set npo_line1 = ""
set npo_line2 = ""
set calculate_subvolume = ""
set mask_subvolume_PDB = ""
set mask_subvolume_PDB_file = ""
set mask_subvolume_PDB_radius = ""
set merge_alsoevenodd = ""
#
#$end_vars
#
echo "<<@progress: 1>>"
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
set imagenumber = 1001
set imagename = "merge"
#
echo imagename = ${imagename}
echo imagenumber = ${imagenumber}
#
set scriptname = 2dx_generateMergeMap
set merge_modus=3D
#
source ${proc_2dx}/2dx_merge_makedirs 
#
\rm -f LOGS/${scriptname}.results
#
set date = `date`
echo date = ${date}
#
set IAQP2 = 0
set IVERBOSE = 1
set phastepnum = 1
set phastep = 0.1
#
# The following is to make sure that for the next "Import Images", the default is correctly initialized.
set initialization_reset = "y"
set initialization_executable = "y"
echo "set initialization_reset = ${initialization_reset}" >> LOGS/${scriptname}.results
echo "set initialization_executable = ${initialization_executable}" >> LOGS/${scriptname}.results
#
if ( ${ALAT} == "0" || ${ALAT} == "0.0" ) then
  ${proc_2dx}/protest "ALAT is not defined."
endif
set ALATnew = `echo ${ALAT} | awk '{ if ( $1 < 0 ) { s = -$1 } else { s = $1 }} END { print s }'`
if ( ${ALAT} != ${ALATnew} ) then
  set ALAT = ${ALATnew}
  echo "set ALAT = ${ALAT}" >> LOGS/${scriptname}.results
endif
#
#???? Check this ???
set SCL = 1
echo "SCL = ${SCL}"
#
# contrast for grey plot
set scale = 1
echo "scale = ${scale}"
#
${proc_2dx}/linblock "Using the most recent merged dataset (from register 0)."
#
#############################################################################
${proc_2dx}/linblock "sourcing sym2spsgrp_sub.com"
#############################################################################
#
source ${proc_2dx}/2dx_sym2spcgrp_sub.com
#
set ABANG = `echo $realang | awk '{s=180-$1} END {print s}'`
echo ABANG = ${ABANG}
#
echo SYM = ${SYM}
echo spcgrp = ${spcgrp}
echo CCP4_SYM = ${CCP4_SYM}
#
echo "<<@progress: +10>>"
#
# contrast for grey plot
set scale = 1
echo "scale = ${scale}"
#
set voldim = 250
set voldim_m1 = `echo ${voldim} | awk '{ s = $1 - 1 } END { print s }'`
set voldim_t142 = `echo ${voldim} | awk '{ s = $1 * 1.42 } END { print s }'`
set voldim_t2 = `echo ${voldim} | awk '{ s = $1 * 2 } END { print s }'`
set voldim_t2m1 = `echo ${voldim} | awk '{ s = $1 * 2 - 1 } END { print s }'`
#
set cellx = `echo ${realcell} | cut -d\, -f1`
set cellxm1 = `echo ${cellx} | awk '{ s = $1 - 1 } END {print s}'`
echo ":cellx = ${cellx},    cellxm1 = ${cellxm1}"
#
set celly = `echo ${realcell} | cut -d\, -f2`
set cellym1 = `echo ${celly} | awk '{ s = $1 - 1 } END {print s}'`
echo ":celly = ${celly},    cellym1 = ${cellym1}"
#
set ALATm1 = `echo ${ALAT} | awk '{ s = $1 - 1 } END {print s}'`
echo ":ALAT = ${ALAT},      ALATm1 = ${ALATm1}"
#
#############################################################################
#############################################################################
#############################################################################
# Here for 3D volume generation:
#############################################################################
#############################################################################
#############################################################################
#
echo "<<@progress: +10>>"
#
if ( ${merge_alsoevenodd} == "y" ) then
  #############################################################################
  ${proc_2dx}/linblock "2dx_processor.exe - to transform HKL file into volume for EVEN images"
  #############################################################################    
  echo "# IMAGE: APH/latfitted_even.hkl <HKL: Input HKL EVEN [H,K,L,A,PHI,FOM]>" >> LOGS/${scriptname}.results
  ${bin_2dx}/2dx_processor.exe --hklin APH/latfitted_even.hkl --mrcout volume_even.map -X ${cellx} -Y ${celly} -Z ${ALAT} --gamma ${realang} --res ${RESMAX} -s ${SYM_NAME}
  echo "# IMAGE: volume_even.map <MAP: final map 3D EVEN images>"  >> LOGS/${scriptname}.results
  #
  #############################################################################
  ${proc_2dx}/linblock "extend the map to 2x2x1 unit cells"
  #############################################################################  
  source ${proc_2dx}/2dx_extend_map.com volume_even.map extended_even.map
  echo "# IMAGE: extended_even.map <MAP: final map 3D, 2x2 unit cells, EVEN>"  >> LOGS/${scriptname}.results
  #
  if ( ${calculate_subvolume}x != "0x" ) then 
    source ${proc_2dx}/2dx_create_subvolume.com ${calculate_subvolume} extended_even.map ${realcell} ${ALAT} 3D_sub_even.map
    echo "# IMAGE-IMPORTANT: 3D_sub_even.map <MAP: 3D sub map, EVEN>" >> LOGS/${scriptname}.results
  endif
  #
  echo "<<@progress: +5>>"
  echo "<<@evaluate>>"
  #############################################################################
  ${proc_2dx}/linblock "2dx_processor.exe - to transform HKL file into volume for ODD  images"
  #############################################################################    
  echo "# IMAGE: APH/latfitted_odd.hkl <HKL: Input HKL ODD  [H,K,L,A,PHI,FOM]>" >> LOGS/${scriptname}.results
  ${bin_2dx}/2dx_processor.exe --hklin APH/latfitted_odd.hkl  --mrcout volume_odd.map  -X ${cellx} -Y ${celly} -Z ${ALAT} --gamma ${realang} --res ${RESMAX} -s ${SYM_NAME}
  echo "# IMAGE: volume_odd.map <MAP: final map 3D ODD images>"  >> LOGS/${scriptname}.results
  echo "<<@progress: +5>>"
  #
  #############################################################################
  ${proc_2dx}/linblock "extend the map to 2x2x1 unit cells"
  #############################################################################  
  source ${proc_2dx}/2dx_extend_map.com volume_odd.map extended_odd.map
  echo "# IMAGE: extended_odd.map <MAP: final map 3D, 2x2 unit cells, ODD>"  >> LOGS/${scriptname}.results
  #
  if ( ${calculate_subvolume}x != "0x" ) then 
    source ${proc_2dx}/2dx_create_subvolume.com ${calculate_subvolume} extended_odd.map ${realcell} ${ALAT} 3D_sub_odd.map
    echo "# IMAGE-IMPORTANT: 3D_sub_odd.map <MAP: 3D sub map, ODD>" >> LOGS/${scriptname}.results
  endif
  #
  #############################################################################
  ${proc_2dx}/linblock "2dx_correlator.exe - for FSC"
  #############################################################################  
  ${bin_2dx}/2dx_correlator.exe --vol1 volume_even.map --vol2 volume_odd.map -R ${RESMAX} --bins 100 --fsc even_odd_fsc.dat
  #
  ${app_python} ${proc_2dx}/plotFSC.py even_odd_fsc.dat PS/even_odd_fsc.ps
  #
  echo "# IMAGE: even_odd_fsc.dat <TXT: FSC data between even odd datasets>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: PS/even_odd_fsc.ps <PS: FSC plot between even odd datasets>" >> LOGS/${scriptname}.results
  #
  echo "<<@progress: +10>>"
  echo "<<@evaluate>>"
endif
#
#############################################################################
${proc_2dx}/linblock "2dx_processor.exe - to transform HKL file into volume for all  images"
#############################################################################    
echo "# IMAGE: APH/latfitted.hkl <HKL: Input HKL all  [H,K,L,A,PHI,FOM]>" >> LOGS/${scriptname}.results
${bin_2dx}/2dx_processor.exe --hklin APH/latfitted.hkl --mrcout volume.map -X ${cellx} -Y ${celly} -Z ${ALAT} --gamma ${realang} --res ${RESMAX} -s ${SYM_NAME}
#
echo "# IMAGE: volume.map <MAP: final map 3D>"  >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "extend the map to 2x2x1 unit cells"
#############################################################################  
echo "<<@progress: +10>>"
source ${proc_2dx}/2dx_extend_map.com volume.map extended.map
echo "# IMAGE: extended.map <MAP: final map 3D, 2x2 unit cells>"  >> LOGS/${scriptname}.results
if ( ${calculate_subvolume}x != "0x" ) then 
    echo "<<@progress: +10>>"
    source ${proc_2dx}/2dx_create_subvolume.com ${calculate_subvolume} extended.map ${realcell} ${ALAT} 3D_sub.map
    echo "# IMAGE-IMPORTANT: 3D_sub.map <MAP: 3D sub map>" >> LOGS/${scriptname}.results
endif
#
#############################################################################
${proc_2dx}/linblock "2dx_processor.exe - to transform HKL file into volume for PSF"
#############################################################################  
#
echo "<<@progress: +10>>"
#
${bin_2dx}/2dx_processor.exe --hklin APH/latfitted.hkl --mrcout PSF.map -X ${cellx} -Y ${celly} -Z ${ALAT} --gamma ${realang} --res ${RESMAX} -s ${SYM_NAME} --psf
#
echo "# IMAGE: PSF.map <MAP: PSF>"  >> LOGS/${scriptname}.results
#
echo "<<@evaluate>>"
#
#
#############################################################################
#############################################################################
#############################################################################
${proc_2dx}/linblock "Now for 2D processing:"
#############################################################################
#############################################################################
#############################################################################
#
#
echo "<<@progress: +10>>"
#
${bin_2dx}/2dx_projector.exe --mrcin extended.map --mrcout ${imagename}-${SYM}.mrc --axis z
#
echo "# IMAGE-IMPORTANT: ${imagename}-${SYM}.mrc <${SYM}-symmetrized final 2D map>" >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "labelh - to transform the map into BYTE format with automatic scaling"
#############################################################################
#
echo "<<@progress: +10>>"
#
${bin_2dx}/labelh.exe << eot
${imagename}-${SYM}.mrc
-3
SCRATCH/${imagename}-${SYM}.mrc
eot
#
#############################################################################
${proc_2dx}/linblock "mrc2tif - to create a TIFF file from the map that can be read by SPIDER"
#############################################################################
#
if ( ! -d RESULTS-TIFF ) then
  \mkdir RESULTS-TIFF
endif
#
${bin_2dx}/mrc2tif.exe << eot
SCRATCH/${imagename}-${SYM}.mrc
RESULTS-TIFF/${imagename}-${SYM}.tif
Y
eot
#
echo "# IMAGE: RESULTS-TIFF/${imagename}-${SYM}.tif"  >> LOGS/${scriptname}.results
${proc_2dx}/linblock "TIFF file created: RESULTS-TIFF/${imagename}-${SYM}.tif"
#
echo "<<@progress: 100>>"
${proc_2dx}/linblock "Done."
#
#############################################################################
#
